The Urbach tail in silica glass from first principles 
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We present density-functional theory calculations of the optical absorption spectra of silica glass for temper- 
atures up to 2400 K. The calculated spectra exhibit exponential tails near the fundamental absorption edge that 
follow the Urbach rule, in good agreement with experiments. We also discuss the accuracy of our results by 
comparing to hybrid exchange correlation functionals. By deriving a simple relationship between the exponen- 
tial tails of the absorption coefficient and the electronic density-of-states, we establish a direct link between the 
photoemission and the absorption spectra near the absorption edge. This relationship is subsequently employed 
to determine the lower bound to the Urbach frequency regime. Most interestingly, in this frequency interval, the 
optical absorption is Poisson distributed with very large statistical fluctuations. Finally, We determine the upper 
bound to the Urbach frequency regime by identifying the frequency at which transition to Poisson distribution 
takes place. 

PACS numbers: 78.40.Pg, 78.20.Bh, 71.23.An, 71.15.Mb 
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At finite temperatures, the absorption spectra of insulators 
can be modified substantially through interaction of the elec- 
tronic states with lattice vibrations. In 1953 Urbach Cll ob- 
served an exponential energy dependence of the absorption 
coefficient near the fundamental absorption edge that varied 
with temperature as follows 



a(w, T) = ao exp 



hu>o(T) — fiui 
kf 



(1) 



Here ljq(T) is a linear function of temperature, which at zero 
Kelvin is defined to be the optical gap, and a and ao are con- 
stants that can be extracted from experiments. The so-called 
Urbach rule described by Eq. ((T), has been observed univer- 
sally in crystals as well as glasses, in both semiconductors 
and insulators. Its origin has been discussed extensively over 
the years J2-01 which has lead to the consensus that it arises 
from transitions between localized electronic levels resulting 
from temporal fluctuations of the band-edge electrons into the 
band gap and extended band-like states. Since Urbach be- 
havior involves electron localization 01 in the presence of 
vibrations, it is not obvious that standard electronic structure 
theories such as the density-functional theory (DFT) in the 
local-density (LDA) or gradient-corrected (GGA) approxima- 
tions can capture these effects, in particular when considering 
their systematic underestimation of band gaps of insulators. 
In a pioneering work, Drabold et al. performed ab-initio LDA 
molecular-dynamics (MD) simulations of amorphous Si and 
calculated the fluctuations in the single-particle energies at 
the band edges [0). These were found to be in good agree- 
ment with the band tail widths deduced from photoemission 
spectra demonstrating the applicability of ab-initio MD in the 
adiabatic approximation to describe the electronic structure of 
semiconductors at finite temperatures. The quantitative com- 
putation of the Urbach tail of the optical absorption, however, 
remains a daunting task. It requires calculating the probabil- 
ity distribution of rare dipole transition events, which neces- 
sitates long-time ab-initio MD simulations in order to obtain 
the time-dependent fluctuations of single-particle energies and 



dipole matrix elements. In the past, almost invariably the as- 
sumption has been made that the latter do not vary appreciably 
with atomic displacements or frequency, an assumption that is 
very difficult to justify especially at high temperatures. 

In this paper, we investigate the Urbach rule in defect-free 
silica glass using ab-initio MD simulations. This has been 
motivated by the need to develop a better understanding of 
the process of laser damage to silica optics, which is of im- 
portance to diverse fields ranging from opto-electronics to in- 
ertial confinement fusion. In the past, a considerable amount 
of computational studies has been directed toward understand- 
ing zero Kelvin absorption due to defects in silica Re- 
cently, the role of temperature has been emphasized by ex- 
periments where damage was generated far below the bulk 
material threshold by photons of energy 3.55 eV when silica 
was heated to about 2200 K 111 111 . The Urbach rule plays a 
crucial role here since the exponential dependence of absorp- 
tion on temperature in Eq. (Q3, necessitates the existence of 
a critical temperature T c , at which the glass absorbs more 
photon energy than it can dissi pate leading to thermal run- 
away and macroscopic damage 111 111 . However, extrapolation 
to higher temperatures of the experimental spectra available 
up to 1900K 11211 predicts a T c that is several hundred de- 
grees higher than the measured value. Therefore better un- 
derstanding of the kinetics of absorption at finite temperatures 
in the Urbach regime is needed. The objective of this Let- 
ter is to study the physical processes that lead to absorption 
in a temperature and energy range for which experiments are 
not available. The key finding is that in the Urbach regime 
absorption occurs locally in space as a Poisson process at sub- 
terahertz frequencies. On the basis of this work, models that 
enable quantitative predictions of T c can be obtained. 

The MD simulations presented in this work are per- 
formed within the DFT-GGA framework using the PW91 
parametrization lfl3i [3 as implemented in the Vienna ab- 
initio simulation package lfl5ll using the projector augmented 
wave method llal . All calculations involve supercells con- 
taining 24 Si02 formula units and the Brillouin zone is sam- 
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FIG. 1. (a) Density-of-states for the perfect glass at zero Kelvin cal- 
culated using the PBE0 hybrid functional (lowermost line) and the 
PW91 functional with scissors correction (PW91+A 9 , second line 
from the bottom). Also shown is the density of states obtained from 
MD simulations at three different temperatures using PW91+A 9 . (b) 
Single-particle energies of the four topmost valence band states and 
the two lowermost conduction band state as a function of time along 
an MD trajectory obtained using two different hybrid functionals 
(HSE06, PBE0) in comparison with the PW91+A 9 approach. 



pled by a 2 x 2 x 2 Monkhorst-Pack fc-point grid. In order 
to obtain a realistic glass model, we started from a liquid 
silica model obtained previously lfl7ll . which was quenched 
down to zero Kelvin over a period of about 10 ps. The ex- 
amination of the electronic structure of the resulting configu- 
rations revealed defect states due to the presence of stretched 
and broken bonds. The defect states were eliminated from 
the model by optimization via a bond-switching Monte Carlo 
(BSMC) technique 0]. Several BSMC-refined configura- 
tions were generated, each representing a random network 
with fixed bond lengths and angles. Subsequently, the con- 
figurations were structurally relaxed to the local GGA total 
energy minimum. The final glass model that was chosen from 
this set had preserved its bond lengths and angles after the re- 
laxation process. The electronic density-of-states (DOS) of 
this configuration is shown in Fig. [TJa) in comparison with 
the PBE0 hybrid functional, which incorporates 25% exact 
exchange 113, 2(J ■ The GGA DOS includes a band gap shift 



of A g = 2.6 eV to account for the systematic underestimation 
of the band gap. The excellent agreement between the two 
calculations is not surprising since in general it is expected 
that GGA-DFT, with the exception of the the band gap, can 



reproduce most of the features of the electronic structure of 
insulators accurately. However, as mentioned earlier, the elec- 
tronic states at the band edges undergo localization when they 
fluctuate inside the band gap. If the magnitude of the GGA- 
DFT band gap shift depends strongly on the degree of local- 
ization of these states, then it would not be possible to develop 
a quantitative prediction of the Urbach tail with this level of 
theory. In order to address this issue, we have explicitly com- 
pared the time evolution of the band edge states at 2200 K, 
calculated from GGA with PBE0 as well as HSE06 (H (a hy- 
brid functional that includes screened exchange) calculations, 
see Fig. 01t>). It appears that with a constant band gap shift 
(2.6 eV for GGA and 0.9 eV for HSE06), all calculations can 
be brought in agreement with each other. We thus conclude 
that DFT-GGA provides a reasonable basis for modeling the 
Urbach tail from first principles. A similar conclusion was 
reached by Alkauskas et al. in studying point defects ifioll . 

The absorption coefficient for photons of energy Tioj of an 
atomic configuration X, can be calculated as follows 



(2) 



where e(ui: X) is the complex dielectric function e = e R+iej. 
In the velocity gauge, ej can be directly computed from the 
single-particle wave functions and energies If22[ f23"ll as well as 
their occupancies as follows 



ej(w;X) = 



4n z e 



2^2 



,2 



x <S(A 9 



(U k -f nk )\M^ n ,(X)\ 
e„'k(X) — e n k{X) — hw) , (3) 



where M^ n ,(X) are the polarization-averaged dipole matrix 
elements between the states nk in the valence band and n'k in 
the conduction band. The sums in Eq. (O run over bands and 
spins, and the real part en can be obtained from ej through a 
Kramers-Kronig relation. Since the latter involves an integra- 
tion over the entire frequency spectrum, we have included as 
many as 1000 unoccupied bands in our calculations in order 
to obtain accurate values for €r. 

At finite temperatures, the response functions as well as 
the DOS are calculated by classical ensemble averaging over 
ionic displacements in the Born-Oppenheimer approximation, 
which amounts to averaging over the MD simulation time 
steps. In this way, the electronic transitions are treated as 
instantaneous. The finite-temperature electronic state occu- 
pancies are determined by the Fermi-Dirac distribution = 
1/ [1 + exp (e n k — /x(T)//csT)]. The electron chemical po- 
tential /i(T) is calculated from the charge neutrality condition, 

J- (/°( £ ))t ^ e = wnere N e is the total number of 
electrons, and (p(e)) T is the average DOS. FigureQJa) shows 
the average DOS at three different temperatures. The gray 
region depicts the band gap narrowing with increasing tem- 
perature, which contributes to creating holes and electrons in 
the valence and the conduction bands. The equilibrium con- 
centrations of free electrons as a function of temperature can 
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be calculated by summing up the total occupancies of the con- 
duction band states. Although the free electron concentration 
can be as large as 10 17 cm -3 at 2400 K, it is still too small 
to have any measurable impact on the absorption coefficients 
in the Urbach regime. The latter can thus be calculated us- 
ing zero Kelvin occupancies and thus neglecting free-electron 
absorption. The DOS is connected to absorption through the 
joint density-of-states (JDOS), which for zero Kelvin occu- 
pancies is defined as J(uj) = J p v (ui')p c (u' + uj) dui', where 
p v ( c \ (lo) is the DOS of the valence (conduction) bands. At fi- 
nite temperatures, a direct relationship between the JDOS and 
the DOS only exists if the fluctuations in the valence and the 
conduction bands are independent 



(p v (uj')} T (p c {u)' + u)) T duj'. 



(4) 



We find that the above is indeed a very good approximation 
over the entire frequency range as illustrated in Fig. EJa), 
which shows the structure of the low-energy exponential tail 
of the JDOS for two temperatures. The figure also illustrates 
that the low-frequency exponential tails of dielectric function 
and JDOS coincide when the latter is scaled by a temperature- 
dependent effective dipole transition probability ~p(T), 



{ej{w)) T « /Z(T) (J( U )), 



(5) 



This important result suggests that exact knowledge of the 
matrix elements is not necessary to reproduce the frequency 
dependence of the low-energy exponential tail of the dielec- 
tric function. They cannot, however, be neglected entirely, 
since they lead to an effective temperature-dependent coeffi- 
cient, the inverse of which decreases linearly with temperature 
and changes by a factor of two between 1200 and 2400 K, as 
shown in the inset of Fig. Efa). 

After obtaining the dielectric functions at different temper- 
atures, it is straightforward to calculate the absorption spectra 
according to Eq. (O, which leads to the curves in Fig. |2jb). 
Fitting the low-energy exponential tails of these spectra to 
Eq. (Q]i, yields a linear temperature dependence for wrj(T) 
in accordance with Urbach's rule, as shown in the inset of 
Fig. EJb). Furthermore, we obtain a = 0.473 in fair agree- 
ment with the experimental value of a = 0.585 reported by 
Saito and Ikushima IU2I1 . The data also enable us to calculate 
the Tauc gap, which is the threshold energy for the onset of 
extended-to-extended electronic transitions. Following Saito 
and Ikushima III 211 who defined the Tauc gap as the photon en- 
ergy corresponding to 5 x 10 3 cm _1 absorption, one obtains 
the data shown in the inset of Fig. |2jb) demonstrating excel- 
lent agreement. The calculated linear temperature dependence 
of ujo(T) in Fig.|2jb) extends the Urbach rule up to 2400 K. 

It is interesting to note that in the Urbach regime the ratio 
(e/(w)) T / (en(uj)) T ^ 1. A first order Taylor expansion 
of Eq. (O with respect to this quantity yields the following 
expressions for the absorption coefficient 



(J(w)) T . (6) 
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FIG. 2. (a) JDOS rescaled by Jl(T) (solid lines) and imaginary di- 
electric function (symbols) at two different temperatures calculated 
from MD simulations. The inset shows ~p(T)~ 1 as a function of tem- 
perature, (b) Absorption coefficient in the Urbach regime at different 
temperatures calculated from MD simulations using the PW91 func- 
tional with a scissors shift. The inset shows the result of the fit to 
Eq. QJ and the comparison with the experimental data from Saito 
and Ikushima fI3l . The shaded region depicts the Urbach regime. 



The second approximation above is obtained by a zeroth order 
expansion about the static dielectric constant (e/?(0)) T , and 
utilizes our earlier finding that the Urbach tail of the imagi- 
nary dielectric function can be obtained from the JDOS, see 
Eq. Eq. (|5). Using Eqs. Eq. © and Eq. ©, we can establish a 
simple relationship between the absorption coefficient and the 
DOS in the vicinity of the absorption edge, where the temper- 
ature dependence of the prefactor mainly enters through the 
effective dipole-transition probability ~p(T), while (en(0)) T 
varies only weakly with temperature, i.e. from 1.81 at OK 
to 1 .99 at 2400 K, corresponding to an increase of less than 
10%. We point out that the DOS is obtained experimentally 
through photoemission spectroscopy. Therefore, the above re- 
sult provides a direct link between the photoemission and the 
optical absorption experiments in the Urbach tail region of the 
spectrum. 

Equation © has the important implication that at finite tem- 
peratures the dipole matrix elements average out such that 
in the Urbach tail they can be replaced by a temperature- 
dependent prefactor. The atomic vibrations also lead to sig- 
nificant statistical fluctuations, which at present are computa- 
tionally too expensive for us to quantify accurately. However, 
we can use Eq. (O to obtain a lower bound on the statistical 
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FIG. 3. (a) First and second moment of the joint density of states 
(JDOS) as well as their ratio at 2400 K. (b) In the Urbach regime 
The standard deviation of the JDOS equals the square root of the 
average indicative of a Poisson distribution. The energy at which 
the ratio deviates from one provides a simple measure for the upper 
limit of the Urbach regime as indicated by the arrow, (c) At low 
energies, the JDOS exhibits two distinct regions corresponding to 
localized-delocalized (Urbach regime) and localized-localized tran- 
sitions. Each region is described by a different exponential the cross- 
ing point of which defines a lower limit for the Urbach regime. Note 
that the Tauc gap lies above the region within which the JDOS ex- 
hibits exponential tails, (d) Upper and lower limits of the Urbach 
regime extracted from the data presented in (b) and (c). 



fluctuations in a(ui) by neglecting the fluctuations in the ma- 
trix elements and study (j7~ 2 (oj)) t . This is shown in Fig.|3a), 
where comparison is made with (^J(oS)) T . A significant de- 
parture between the two curves is observed in the Urbach tail, 
where the ratio ( J' 2 (uj)) t / (J'(uj)) t can reach 10 4 . Hence 
the standard deviation from mean absorption for frequencies 
in the Urbach tail is up to two orders of magnitude larger than 
(a(uj)) T . 

The large fluctuations in the Urbach regime originate from 
the discrete nature of the JDOS itself. Even if (J(u>)) T <C 1, 



there can only exist an integer number Af of pairs of electronic 
states available for transition at each instant of time. Hence 
the absorption coefficient locally fluctuates between zero and 
/I 2 x Af, which can amount to fluctuations much larger than 
the mean value (J) T itself. This behavior is a consequence 
of the quantum nature of matter. Let us define an absorption 
event as the instant of time when Af > 0. Whenever absorp- 
tion events occur so rarely that they can be considered inde- 
pendent, the absorption process is Poisson distributed, with 
(J(uj)) t interpreted as the average rate of occurrence. An 
important signature of the Poisson distribution is that its stan- 
dard deviation is equal to the square root of its average. As 
shown in Fig.[3jb) in the Urbach regime, this equality indeed 
holds while the ratio rapidly decreases at only slightly higher 
energies. The sharp transition provides a natural definition for 
the frequency uju which signifies the upper bound to the Ur- 
bach region and below which the statistics become Poisson 
distributed. Figure [3jd) shows that in the temperature range 
considered here, uijj decreases linearly. These observations 
bestow the finite-temperature JDOS in the Urbach tail with 
distinct physical significance as it represents the average rate 
of occurrence of absorption events in this frequency interval. 

The Urbach tail region is a closed frequency interval in 
the optical spectrum. Above, we have determined the upper 
bound for this interval. We can also determine its lower bound 
using Eq. (HJi to compute (J(w)) T down to very small values 
with good statistical accuracy. Figure [3c) shows that at very 
low energies there is a transition in the signature exponen- 
tial decay of the Urbach tail to a steeper decay curve. The 
frequency at which this transition occurs can be used to 
identify the lower bound of the Urbach region, the tempera- 
ture dependence of which is shown in Fig. Od). The optical 
transitions in this lower-frequency region correspond to tran- 
sitions between localized levels in the exponential tails of the 
valence and the conduction bands while the Urbach tail orig- 
inates from transitions between localized tail states and ex- 
tended band-like states. 

Electronic localization in the Urbach tail has been discussed 
extensively in the literature As will be described in 

detail in a future publication, we indeed observe pronounced 
localization in our simulations. Here it suffices to point out 
that in the Urbach region the conduction band edge is prac- 
tically flat in reciprocal space which is is a manifestation of 
the localization of these eigenstates within the 72-atom super- 
cell. In contrast, the conduction band minimum of the perfect 
glass exhibits significant dispersion in reciprocal space with a 
standard deviation of 250 meV. 

Finally, let us discuss the impact of the above findings on 
the modeling of laser heating in silica. Commonly, this pro- 
cess is described by a heat conduction equation with a source 
term that incorporates energy deposition by linear coupling to 
the laser light, a ((J, T) I(v), where /(r) is the laser light in- 
tensity. Neglecting fluctuations, this term can be parametrized 
by T) = (a(w)) T . However, the rare event nature of ab- 
sorption in the Urbach regime calls for a(u>, T) to be treated 
as a discrete Poisson process, where at a rate proportional to 
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(J(u})) T an absorption event of strength iJJi 2 jcU (ejt(0)) T 
takes place. The large fluctuations introduced in this way can 
reduce the predicted thermal run-away temperature leading to 
better agreement with experiments. 

This work performed under the auspices of the U.S. Depart- 
ment of Energy by Lawrence Livermore National Laboratory 
under Contract DE-AC52-07NA27344 with support from the 
Laboratory Directed Research and Development Program. 
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